function [evn, isamp, icoff, crit, template, predelay] = ClementsBekkers(data, samplerate, rise, decay, threshold, sign, dispflag, lpfilter, template_type, ntau, hmini)
%
% implement Clements and Bekkers algorithm
% Biophysical Journal, 73: 220-229, 1997.
%
% incoming data is raw data (or filtered)
% rise and decay determine the template that 
% is the matching template to work against the data at each point
% threshold sets the detection limit.
% the returned eventlist is a list of indices to the events detected.
%
% This code uses the optimizations discussed in the article.
%
% 19 June 2002 Paul B. Manis, Ph.d. pmanis@med.unc.edu
%%


evn = [];
isamp = [];
icoff = [];
crit = [];
switch(nargin)
    case 0 % testing mode....
        [data, samplerate, rise, decay, threshold, sign, dispflag] = testdata;
        fsamp = 1000/samplerate; % get sampling frequency
        fco = 1600;		% cutoff frequency in Hz
        wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
        if(wco < 1) % if wco is > 1 then this is not a filter!
            [b, a] = butter(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
            data = filter(b, a, data); % filter all the traces...
        end
        sign = 0; % negative detection
        dispflag = 1;
        lpfilter = 1600;
        template_type = 1;
        ntau = 5;
        hmini = [];
    case {1,2,3,4,5}
        QueMessage('ClementsBekkers - insufficient arguments');
        return;
    case 6
        dispflag = 0;
        lpfilter = 1;
        template_type = 1;
        ntau = 5;
        hmini = [];
    case 7
        lpfilter = 1;
        template_type = 1;
        ntau = 5;
        hmini = [];
    case 8
        template_type = 1;
        hmini = [];
    case 9
         hmini = [];

    otherwise
end;
%fprintf(1, 'ntau = %f\n', ntau);

sumy=0;
sumy2=0;
evn=[];
crit=zeros(length(data),1);
cx=zeros(length(data),1);
scale=zeros(length(data),1);
pkl=zeros(length(data), 1);
eventlist=zeros(length(data), 1);
nevent = 0;
[template, predelay] = cb_template(template_type, samplerate, rise, decay, lpfilter, ntau);
N = length(template);
if(template_type == 4) % data
    Npost = length(template);
else
    Npost = floor(decay*ntau/samplerate);
end;
%%
[pkl, eventlist, crit, scale, cx, nevent] = clembek(template, data, threshold, sign);
%fprintf(1, 'clembek says: %d events\n', nevent);
%sume = sum(template); % only need to compute once.
%sume2 = sum(template.^2);
x=find(eventlist > 0);
nevent = length(x);
eventlist = eventlist(x);
pkl = pkl(x);
% Measurement is done in C code in clembek. The following is just matlab
% origina used to develop the C code. 
% for i = 1:length(data)-N
%   if(i == 1)
%      sumy = sum(data(1:N)); % optimization...
%      sumy2 = sum(data(1:N).^2);
%   else
%      sumy = sumy + data(i+N-1) - data(i-1);
%      sumy2 = sumy2 + (data(i+N-1)^2) - (data(i-1)^2);
%   end;
%   sumey = sum(data(i:N+i-1).*template);
%   S = (sumey - sume * sumy/N)/(sume2 - sume * sume/N);
%   C = (sumy - S*sume)/N;
%   f = S*template+C;
%   SSE = sumy2 + (S*S*sume2)+(N*C*C)-2*(S*sumey + C*sumy-(S*C*sume));
%   CRITERIA = S/sqrt(SSE/(N-1));
%  crit(i) = CRITERIA;
%  scale(i) = S;
%  cx(i) = C;
%  if(sign == 0 & CRITERIA <= -threshold)
%      eventlist = [eventlist i];
%      [mpk ipk] = min(data(i:N+i-1));
%      pkl = [pkl ipk+i-1];
%   end;
%    if(sign == 1 & CRITERIA >= threshold)
%      eventlist = [eventlist i];
%      [mpk ipk] = max(data(i:N+i-1));
%      pkl = [pkl ipk+i-1];
%   end;

%end;
% now reduce the eventlist to single points. 
if(nevent == 0)
    QueMessage('ClementsBekkers: No Events Detected', 1);
    return;
end;

nev = 0;
evl = []; % NaN*zeros(1, nevent);
evn = evl;
pks = evl;
isamp = evl;
icoff = evl;

%%
for i = 1:nevent
    if(i < nevent && (eventlist(i+1)-eventlist(i)) == 1)
        evl = [evl eventlist(i)]; % build current event indices in evl array
    elseif(i == nevent)
        evl = [evl eventlist(i)];

        if(~isempty(evl))
            nev = nev + 1;
            if(sign == 0)
                [emin, imin] = min(crit(evl)); % end of current event - find best crit/fit
            else
                [emin, imin] = max(crit(evl)); % end of current event - find best crit/fit
            end;
            
            ipt = imin+evl(1)-1;
            %      fprintf('event %d: imin = %f emin = %f crit = %f amp = %f\n', nev, ipt, emin, crit(ipt), s(ipt));
            evn = [evn ipt]; % store index here...
            pks = [pks pkl(i)];
            isamp = [isamp scale(ipt)];
            icoff = [icoff cx(ipt)];
        end;
        evl = []; % clear array for next event.
    else
        if(~isempty(evl))
            nev = nev + 1;
            if(sign == 0)
                [emin, imin] = min(crit(evl)); % end of current event - find best crit/fit
            else
                [emin, imin] = max(crit(evl)); % end of current event - find best crit/fit
            end;
            ipt = imin+evl(1)-1;
            %       fprintf('event %d: imin = %f emin = %f crit = %f amp = %f\n', nev, ipt, emin, crit(ipt), s(ipt));
            evn = [evn ipt]; % store index here...
            pks = [pks pkl(i)];
            isamp = [isamp scale(ipt)];
            icoff = [icoff cx(ipt)];
        end;
        evl = []; % clear array for next event.
    end;
end;
evn = evn + floor(predelay/samplerate)-1; % replace list with new list of best points.
eventlist = eventlist + floor(predelay/samplerate);   
fprintf(1, 'clembek identified %d events\n', length(evn));
%pks = pks + floor(predelay/samplerate)-1;
if(dispflag > 0)
    hf = newfigure('EPSC_minitraces');
    if(dispflag == 1)
        clf;
    end;
    t = samplerate*(0:(length(data)-1));
    h1=subplot('position', [0.07, 0.55, 0.8, 0.4]);
    set(h1, 'Tag', 'cbtestfig1');
    if(dispflag == 2)
        hold on;
    end;
    plot(t, data, 'k');
    hold on;
    %   plot(t(eventlist), data(eventlist), 'rx');
    plot(t(pks), data(pks), 'ro', 'markersize', 4, 'markerfacecolor', 'red');
    plot(t(evn), data(evn), 'go', 'markersize', 4, 'markerfacecolor', 'green');
    set(gca, 'YLim', [min(data) max(data)]);
    h2=subplot('position', [0.07, 0.3, 0.8, 0.2]);
    set(h2, 'Tag', 'cbtestfig2');
    if(dispflag == 2)
        hold on;
    end;
    plot(t, crit);
    hold on
    plot([t(1),t(end)],[threshold,threshold], 'k-');
    plot([t(1),t(end)],-[threshold,threshold], 'k-');
    h3=subplot('position', [0.07, 0.07, 0.8, 0.2]);
    set(h3, 'Tag', 'cbtestfig1');
    if(dispflag == 2)
        hold on;
    end;
    plot(t, scale, 'k');
    hold on
    plot(t, cx, 'r');
    u = get(gca, 'Xlim');
    set(gca, 'Xlim', [0 u(2)]);
end;

return;
%% 

function  [data, samplerate, rise, decay, threshold, sign, dispflag] = testdata()

samplerate = 0.1;
rise = .05;
decay = .3;
sign = 0;
threshold = 4;
dispflag = 1;
data = zeros(1, 12000);
noise = randn(1, length(data));
% filter the noise then put in the data array
data = data + noise;
fsamp = 1000/samplerate; % get sampling frequency
fco = 2000;		% cutoff frequency in Hz
wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
if(wco < 1) % if wco is > 1 then this is not a filter!
    [b, a] = butter(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
    data = filter(b, a, data); % filter all the traces...
end
%noise2 = 6*randn(length(data), 1);
%fco = 30;		% cutoff frequency in Hz
%wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
%if(wco < 1) % if wco is > 1 then this is not a filter!
%   [b, a] = butter(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
%   noise2 = filter(b, a, noise2); % filter all the traces...
%end
%data = data + noise2; % adding low frequency noise
predelay = decay*1;
N = floor((decay*5+predelay)/samplerate);
for i = 1:N
    t = i*samplerate;
    if(t >= predelay)
        template(i) = (1-exp(-(t-predelay)/rise*0.5))*exp(-(t-predelay)/decay);
    else
        template(i) = 0;
    end;
end;
m=max(template);
template = -template/m; % normalize to size 1.

for i = 1:10
    amp = i;
    pos = floor(i*10000/11);
    data(pos:pos+N-1) = data(pos:pos+N-1)+template*i;
end;
return;
